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Abstract. We present an algorithm for two-dimensional 
radiative transfer in axisymmetric, circumstellar media. 
The formal integration of the transfer equation is per- 
formed by a generalization of the short characteris- 
tics (SC) method to spherical coordinates. Accelerated 
Lambda Iteration (ALI) and Ng's algorithm are used to 
converge towards a solution. By taking a logarithmically 
spaced radial coordinate grid, the method has the natural 
capability of treating problems that span several decades 
in radius, in the most extreme case from the stellar ra- 
dius up to parsec scale. Flux conservation is guaranteed 
in spherical coordinates by a particular choice of dis- 
crete photon directions and a special treatment of nearly- 
radially outward propagating radiation. The algorithm 
works well from zero up to very high optical depth, and 
can be used for a wide variety of transfer problems, includ- 
ing non-LTE line formation, dust continuum transfer and 
high temperature processes such as compton scattering. In 
this paper we focus on multiple scattering off dust grains 
and on non-LTE transfer in molecular and atomic lines. 
Line transfer is treated according to an ALI scheme for 
multi-level atoms/molecules, and includes both random 
and systematic velocity fields. The algorithms are imple- 
mented in a multi-purpose user-friendly radiative transfer 
program named RADICAL. We present two example com- 
putations: one of dust scattering in the Egg Nebula, and 
one of non-LTE line formation in rotational transitions of 
HCO + in a flattened protostellar collapsing cloud. 
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1. Introduction 

Molecular line and dust continuum observations are an im- 
portant tool for studying the envelopes and disks around 
young stellar objects (YSO), post-AGB stars and AGN. 
One of the main difficulties in interpreting such obser- 
vations is that optical depths effects play an important 



role in the emission of this radiation. It is, for instance, 
well known that self-absorption and non-LTE effects are 
the main processes at work in shaping the characteristic 
asymmetric double-peaked emiss ion li ne profiles from col- 
lapsing protostellar cores (Zhou, |l992|) . Radiative transfer 
computations, at an appropriate level of complexity, are 
therefore needed in order to reconstruct the density, ve- 
locity and temperature structure of an observed cloud. 

If densities drop below the critical density, the sys- 
tem may deviate from local thermodynamic equilibrium 
(LTE). Line trapping and photon escape from the line 
wings can in some cases be treated in the large-velocity- 
gradient (LVG) limit, but this approach is only valid if 
systematic velocity fields are much greater than the local 
turbulent line width. In all other cases one must perform 
a full non-LTE line transfer computation. 

For problems that can be formulated in 1-D slab or 
spherical geometry there exist many radiative transfer 
programs, many of which use sophisticated techniques 
such as Ac celerat ed Lambda Iteration (ALI; see a review 
by Hubeny 1989 ) and Complete Linearization (CL; Auer 
& Mihalas 1969 ). But the requirement of 1-D geometry is 
often too restrictive. Distinctly non-spherical features are 
often observed fro m YS O, such as bipolar reflection neb- 
ulae (e.g. Lenzen 1987), bipolar outflows (see Bachillcr 



1996[ ) and disks (e.g. McCaughrean & O'Dell |l996|) . Even 
the progenitors of these YSOs, starless dense cloud cores, 
seem to appear as elongated structures at millimeter wave- 
lengths (e.g. Myers et al. 1991 ), indicating that even in the 
early stages of star formation spherical symmetry does not 
apply. The case of post-AGB stars and Planetary Nebulae 
is just as compelling, with the majority of these nebulae 
being bipolar. The Cygnus Egg Nebula (CRL 2688) and 
te Red Rectangle (HD 44179) are perhaps the most spec- 
tacular examples of such bipolarity. 

To model such objects, clearly one must resort to 
multi-dimensional transfer computations. There is a vast 
literature on this topic. Methods roughly fall in one of 
three catagories: Monte Carlo methods, Discrete Ordinate 
methods and Angular Moment methods. Monte Carlo 
codes are very flexible and can be used for a large va- 
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riety of problems in multidimensional geometries, such as 
UV continuum transfer (e.g. Spaans 1996), optical and 
infrared continuum transfer (Wolf et al. 1999), molecu- 
lar line transfer (Hogerheijde 1998] ), and Compton scat- 
tering (e.g. Pozdnyakov, Sobol & Sunyaev 1977; Haardt 
& Maraschi 1991). Such methods perform well at low to 
medium optical depths, but it is well known that at high 
optical depths they converge very slowly. 

Angular Moment methods, on the other hand, are very 
well suited to treat the high optical depths regime, since 
they are related to (or variants of) the diffusion equation 
(see e.g. Yorke et al. 1993, Sonnhalter et al. 1995, Murray 
et al. 1994). However, it is not surprising that they fail at 
low optical depth, since the diffusion approximation was 
never meant for this regime. 

In the Discrete Ordinate approach, not only space is 
discretized into cells, but also the photon propagation di- 
rection. Most multi-dimensional implementations of the 
Discrete Ordinate methods are based on the "Lambda 
Iteration" scheme (e.g. Collison & Fix |l99l|, Efstathiou 



& Rowan-Robinson 1991, Philips 1999). The advantage 
of these methods over the Monte Carlo approach is that 
they do not involve random noise, and therefore provide 
'cleaner' answers. But they suffer from the same conver- 
gence problems as Monte Carlo methods. However, for 
Lambda Iteration there are various ways to cure this dis- 
ease. The most well known of these methods is Acceler- 



ated Lambda Iteration (ALI, e.g. Scharmer 1981, Rybicki 
& Hummer pMl]). 



In this paper we will focus on the Discrete Ordinate ap- 
proach to radiative transfer because of its versatility, accu- 
racy and the wide range of convergence acceleration tech- 
niques available. However, despite the relative efficiency 
of these methods, multidimensional calculations remain 
costly. Feasibility constraints can pose severe limits on the 
spatial and angular resolution, which could easily result 
in unacceptable numerical diffusion. Also, this limits the 
number of models one can reasonably make to fit obser- 
vations, which could lead to dangerous undersampling of 
the parameter space. 

The bottleneck lies in the integration of the formal 
transfer equation. The most straightforward way of per- 
forming these integrals is by the method of "Long Char- 
acteristics" , which is accurate, but rather costly in CPU 
time. A more efficient algorithm for doing this in two di- 
mensions is the method of "Short Characterises" (SC; 
Mihalas et al. h978L Kunasz & Auer |1988l, Auer & Pale- 



tou 1994, Stone et al. (1992). These algorithms are de- 
signed specifically with cartesian or cylindrical coordinates 
in mind, and are not straightforward to generalize to other 
coordinate systems. For circumstellar envelopes, however, 
there are several arguments favoring the use of spherical 
(polar) coordinates, as opposed to cylindrical coordinates. 
Most circumstellar nebulae have density and temperature 
profiles that are peaked towards the center. This means 
that the radiation field is dominated by photons emitted in 



the central regions, which are subsequently reprocessed in 
the outer parts of the nebula. The numerical scheme must 
therefore be able to resolve both the very concentrated 
central regions and the extended outer regions simultane- 
ously. Also, it must guarantee that all radiation emitted 
at small radii will eventually emerge at large radii, which 
amounts to saying that flux must be conserved over a large 
range of radii. Using spherical coordinates and a logarith- 
mic radial grid is the most natural way to cover such a 
large dynamic range and guarantee flux conservation. 

The goal of this paper is to describe, test and demon- 
strate an algorithm that generalizes the Short Character- 
istics method to spherical coordinates^] It is an algorithm 
specifically suited for axisymmetric circumstellar nebulae 
and disks. It has been implemented in a multi-purpose ra- 
diative transfer code named RADICAL, which is designed 
to perform 2-D computations in dust continuum emis- 
sion/absorption, multiple scattering off dust grains, non- 
LTE line transfer, and (for application to X-ray binaries 
and AGN) Comptonization. In this paper we describe the 
method of short characteristics in spherical coordinates, 
and focus our attention to the cases of simple isotropic 
scattering off dust grains and non-LTE line transfer for 
multi-level molecules. For a more extensive discussion of 



the algorithm and its applications, sec Dullemond ( 1999| 
The structure of this paper is as follows. In Section 
we will present the equations of transfer we wish to solve. 
In Section^ we shall review the method of short character- 
istics as it is often presented in the literature. In Section 
^ we will show how this method can be generalized to 
spherical coordinates. Then we will put the algorithm to 
the test in Section |^. Finally we will present two example 
applications in Sections ^| and [?]. 

2. Equations of radiative transfer 

The method that will be described in this paper is an 
Accelerated Lambda Iteration method. In such an algo- 
rithm the integration of the formal transfer equation is 
performed using a "Lambda Operator" . In this section we 
will present the equations that are to be solved, and we 
define the Lambda Operator. The numerical details of the 
Lambda Operator will be given in the next sections. 
The formal transfer equation is 



— - = a u (S u - I v ) 
as 



(1) 



with I v is the intensity, S v the source function, a v the 
opacity, and s the path length. This equation must hold 
along every straight line through the medium. Its integral 
form along a ray through a point P reads: 



UP) 



'i AO) 



e-<S{r' v )dr' v 



(2) 



(i 



1 Just prior to subm ission we became aware of a paper by 
Busche & Hillier ( 2000 ) , who describe a method of short char- 



acteristics in spherical coordinates that is quite similar to ours. 
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where r„ is the optical depth along the ray, between point 
P and the edge of the medium. After evaluating this inte- 
gral for all angles u>, one can compute the mean intensity 

Jv 



J" = J~ J Iv{u)du) . 



(3) 



The entire operation of computing J V {P) at every point 
P, for a given source function S v , can be written as the 
action of a linear Lambda Operator A: 



J v = A[S V ] 



(4) 



Using this Lambda Operator we can write down the com- 
plete transfer equation for a simple problem of thermal 
emission and isotropic (dust) scattering 



S v = eB u (T) + (l-e)A[S„], 



where e = a„/a v is the thermalization coefficient (with 
a^ bs the thermal absorption opacity), and B V (T) is the 
Planck function. Solving the transfer problem for isotropic 
scattering and thermal emission amounts to solving Eq. (JbJ) 
for S v . The Lambda Iteration procedure amounts to it- 
eratively applying the Lambda Operator and computing 
the new S v until convergence is reached. The Accelerated 
Lambda Iteration procedure, which converges much faster, 
is a variant of this procedure, involving an approximate 
operator A*. For details we refer to Hubeny (1989) and 



Rutten ( |1999D - 

For multi-level line transfer, we follow the treatment of 
Rybicki & Hummer (1991). Consider an atom or molecule 



having N levels, with spontaneous radiative downward 
transition rates Ay , Einstein coefficients Bij and collision 
rates Cy between levels i and j. The line profile function 
ifij (y) determines at which frequencies the line emits and 
absorbs. When no systematic fluid velocities are present, 
the line profile function is isotropic, and is normalized to 
unity. For the application to circumstellar envelopes, the 
dominant broadening mechanisms are turbulent and ther- 
mal broadening. These two mechanisms produce a Gaus- 
sian profile: 



<Pij (v) = 



O-tot v ij 1 



: exp 



Here c is the speed of light, the line-center frequency 
of the transition between levels i and j, and atot is the line 
width, 



«tot — «turb + 



(7) 



where Tki n is the (kinetic) temperature of the gas, m mo i 
the mass of the molecule, and aturb is the turbulent line 
width. A systematic fluid velocity can cause the line profile 
function to be angle-dependent in the lab frame as a result 
of Doppler shift, 



<Pij(u, v) = <pij(y(l - u> ■ v/c) - Vij) 



(8) 



The opacity in the line associated with this line profile is: 

otij{u),v) = ^-N{njBji -mBijjipijfajv) , (9) 

where rii are the fractional level populations, and N the 
number density of molecules. We assume complete redis- 
tribution for the lines. The source function is then inde- 
pendent of frequency and angle: 



i Ai 



fijBji 



TLiBij 



(10) 



The transfer equation for this source function is then 
dl i3 {u3, v 



ds 



(11) 



(5) where we assume non-overlapping lines. 



The source term Sij is known once the fractional level 
populations rii are known. They are a solution of the sta- 
tistical equilibrium equation. Using the definition of the 
line- integrated Lambda Operator Ay[jSy] = Jy, with 

Jij = J I(u>, u)<pij((jj, v)du)dv , (12) 
the statistical equilibrium equations become: 

^[rijAji + {rijBji - niBij)Aji[Sji]] 
j>i 

~ Yl \- niAl i + ( niBl i ~ "■ l} ^ X h [ S u\] (13) 

3<l 
j 

The non-locality of radiative transfer is now hidden in the 
Aji operator, so that Eq.([l3|) now represents the complete 
(non-linear) set of equations for line transfer. Lambda It- 
eration now proceeds by iteratively applying the Aji op- 
erator and solving the matrix equation represented by 
Eq.(fl3l). Accelerated Lambda Iteration proceeds accord- 
ing to the MALI scheme of Rybicki & Hummer ( 1991 ). 



(6) 

3. Short characteristics in cartesian coordinates 



To carry out the Lambda Iteration or Accelerated Lambda 
Iteration procedure, we need a numerical implementation 
of the Lambda Operator A. In Cartesian coordinates, the 
formal transfer equation Eq.([TJ) becomes 



ds 



dx 



—r- = Ux-z h Uy-— = a v {S v 



(14) 



where translational symmetry in the ^-direction was as- 
sumed. 

The numerical implementation of the Lambda Op- 
erator amounts to integrating Eq.(p^) for given S u and 
a u . This must be done on a 2-dimensional spatial grid 
x = (xi,yj), for a discrete set of directions uj — {uji} and 
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Fig. 1. An illustration of the Long Characteristics (LC) 
method in Cartesian coordinates. The intensity at point 
P is computed integrating the transfer equation along the 
entire ray from the upstream boundary (point U) towards 
point P. 



Fig. 2. An illustration of the SC method in Cartesian co- 
ordinates. The short characteristic is the line connecting 
point U to D through P. The value of the intensity at U is 
determined by quadratic interpolation between the points 
A B and C. 



frequencies v = \v{\. This will provide the specific inten- 
sity I(xi, yj] u>k, w) for all i,j, k, I. Let us focus on a given 
point P = (xi,yj) and on a single direction and frequency, 
u = LJk, v = v\. The integral of Eq.(p^) can be performed 
numerically along the entire characteristic starting at the 
upstream boundary, heading in the downstream direction 
(i.e. the direction where the radiation comes from) and 
ending at point P (see Fig. [l]). This direct approach is 
called the method of Long Characteristics (LC). Provided 
the discretization in angle u> is appropriate, this method 
is quite accurate and reliable. But it has a computational 
redundancy, and hence it is overly time-consuming. Con- 
sider, for instance, a spatial grid N x N, a set of N w di- 
rections and of N v frequencies. The long characteristics 
integral of Eq. ( |l4| ) typically requires in the order of N in- 
tegration steps. This means that while the dimension of 
the grid is N 2 x N u x N v , the total computational time 
scales as 



t C pu oz N x x N v . 



(15) 



The Short Characteristics method of integration (SC; 
see Fig. ||) does not have this redundancy. Instead of per- 
forming the integral along the entire ray (the long charac- 
teristic) , we perform the integral only along that portion of 
the ray (the short characteristic) which connects a point 
U on the grid upstream of P to the closest intersection 
downstream of P itself. The intensity at P is given by 



e-^I v (U-uj) + 



(16) 



where r„ is the optical depth between points U and P. 
The upstream intensity I v (U;uj) can be found from the 
intensities at A, B and C by 3-point quadratic interpola- 
tion 

/„([/; w) = aI v {A;u)+bI v (B;u)+cI v {C;u) (17) 



where a, b and c are the usual Lagrange coefficients for 
polynomial interpolation. Quadratic or higher order inter- 
polation is necessary in order to reproduce the diffusion 
limit for high optical depth, which is governed by a second 
order partial differential equation. 

The integral from U to P can be computed with sec- 
ond order accuracy by interpolating the source function 
S v (x(t{,); ui) between the points D, P and U. Following 
Olson & Kunasz (1987), one finds 



e- T »I„(U;uj)+u„S l/ (U;Lj) 

+p v S u (P; u) + d v S v {D] u>) , 



with 

u v = 

Vv = 

d v = 

eo = 

ei = 

&i = 



eo + [ea - (2t v + f v )ei ]/[t u (t v + f v )] 

\{t v + r v )ei - e-z]/[T v f v ] 

[e 2 - T v e\]/[f v {T v + ?u)} 

1 - e- T " 

t v - e 

tI - 2e, 



(18) 

(19) 
(20) 
(21) 
(22) 
(23) 
(24) 



where t„ and f„ are the depths at U and D, respectively. 
It should be noted that this quadrature formula may have 
pathological behaviour if the source function and/or the 
opacity varies strongly between the points U, P and D. 
This problem can be solved by limiting the resulting in- 
tegrals between zero and max (^(P), j v (JJ)) As, where 
As is the path length along the short characteristic, and 
j u = a v S v . 

By systematically performing the integrals over all the 
short characteristics, one can find an approximate formal 



solution of the transfer equation (Kunasz & Auer 1988 
Auer & Paletou [1994 Auer et al. |1994[ ) . A key ingredient 



for the SC method to work is that the integrals should be 
performed in the right order, so that the upstream inten- 
sities I u {A;<jS), I„(B\u)) and I„(C;u)) are known before 
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Fig. 3. The global and local coordinate systems used to 
describe the radiation field. Here ir/2 — 9 is shown instead 
of 9 for the clarity of illustration. 

the integral is performed and Eq.([l8|) evaluated. In order 
to do so, the grid must be swept from the two upstream 
boundaries towards the two downstream boundaries. 

The method of Short Characteristics is computation- 
ally less time consuming than the method of Long Char- 
acteristics, because now the transfer integral is performed 
over a much shorter path. For the same discretization in- 
troduced earlier in this section, the computational time 
scales as 



t CPU cc N x N u x N v 



(25) 



which is typically an factor N shorter than in the case of 
Long Characteristics. 

4. Short characteristics in spherical coordinates 

We now wish to generalize the Short Characteristics to 
spherical coordinates. In the following we refer to a stan- 
dard spherical coordinate system (R, 0, $) where is the 
latitude and $ the azimuth. By assuming axial symmetry, 
any dependence on <!> is suppressed, although radiation is 
still allowed to travel along 9/9$, as well as in the radial 
and meridional directions. 

In order to describe the radiation field at each spatial 
point P = (R, 0) we need to set up a local coordinate 
system to characterize the photon direction at P. We in- 
troduce two independent angles on the sky of the local 
observer: 9 and <f>. The north pole of this local sky-map 
is chosen to coincide with the outward-pointing radial di- 
rection. The (j> angle is gauged in such a way that (p = 
points parallel to the equator of the global coordinate sys- 
tem (see Fig. U). As is customary in transfer theory, we 
use n = cos 9 instead of 9 itself, so the specific intensity 



depends upon the two spatial variables R, 0, the ray di- 
rection /i, cj), and the frequency v, I = I„(R, ; fi,<f>). The 
transfer equation, Eq.(|l|), in spherical coordinates reads 



dl v dl v . dl v 1 - /i 2 dl v 

ds ~^dR R Sm0 99 + R dn 



i<f> y/1 — (I 2 dl v 

7hf 



(26) 



tan 6 R 



a v (S v 



An important consequence of the use of spherical co- 
ordinates is that, contrary to what happens for cartesian 
coordinates, the photon angles fi and <p are no longer con- 
stant along the rays. The variation of i?, O, \i and <fi along 
the path are 



dR 



fi- 



ds 

dfi 1 — /J 2 
ds 



R 



■ sm < 



dO 
ds 

dcj) \fl — fjp cos ( 



(27) 
(28) 



R ' ds R tan 6 ' 

where s is the path length. Solving these equations yields 

(29) 



R 2 = b 2 + s 2 



cos 6 



sm < 



z + s cos 6 C 

Vb 2 + s 2 ' 
s 



VWTs 2 



b 2 COS Ooo — ZqS 



by/b 2 + s 2 -(z + s cos Goo) 2 ' 



(30) 
(31) 
(32) 



where b is the impact parameter of the ray with respect to 
the origin, z$ is the height above the midplane of closest 
approach to the symmetry-axis, and 0o is the inclination 
at infinity. When projected into the subspace spanned by 
R, the trajectory becomes a hyperbola, as is shown in 
Fig. ^. We stress that this shape is caused by eliminating 
the dependence on the ^-angle, and is purely a projection 
effect. 

For the numerical implementation of the short charac- 
teristics scheme, we are interested in those characteristics 
that pass through a grid point P = (i?fc, 0;) and are tan- 
gent to one of the local discrete ordinates (/ij, 4>j). Clearly, 
once 0;, m, 4>j are fixed, such a characteristic is unique 
and the values of its parameters are 



b 2 = Rl(l-t4 



cos 0oo = (Uj cos 0; + \ 1 — fif sin 0; sin 



zo = i2fc[(l - A*f) COS©; 



(33) 
(34) 



(35) 



The short characteristic passing through (Rk,Qi ; Hi,4>j) 
is defined as the section of this curve that starts at the clos- 
est intersection with the grid lines upstream of P (point 
17), passes through P and ends at the closest intersection 
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Fig. 4. An example of a long characteristic in an axially 
symmetric space. Only the upper quadrant is shown. The 
vertical axis is the symmetry axis and the horizontal axis 
the equator. The bot-dashed lines represent the asymp- 
totes of the hyperbolic characteristic. 

with the grid lines downstream of P (point D). The loca- 
tion of the points U and D is specified by the correspond- 
ing values of parameter s along the ray, sjj and sp>, which 
are found solving equations (^9|) (|3^) with R = Rx and 
6 = Q L , where K = k-l,k, k+1 and L = 1 -1,1,1 + 1. 
Both R = Rk and 9 = 9; need to be included because 
the characteristic may intersect the same 9 or R grid line 
twice. In principle, each equation has two solutions for a 
given value of K and L, giving 12 possible roots 

Sl ... 6 = ±^R K - 6 2 (36) 

S7...12 = — TEr\ -zocos9oo± 

COS 2 9oo — cos 2 <Sl I 

cos9 L y / 6 2 (cos 2 9 00 - cos 2 9 L ) + z 2 J . (37) 

However, two of these solutions always give s = sp, 
i.e. P = (Rk,@i) itself, and are of no interest. Of the 
remaining 10 roots, some are complex and must be re- 
jected. Between the real solutions, the one representing 
point D (U) is selected asking that s > sp (s < sp) and 
that \s— sp \ is minimum. For convenience, in the following 
we will denote with Rx, Oi, fli and 4>i the values of the 
independent variables along the ray at point U. 

Although, as we have just shown, short characteristics 
can be easily defined in spherical coordinates, two major 
problems have to be solved before they can be of any use 
in building a transfer algorithm. The first point concerns 
the fact that, as it was mentioned earlier, fi and <j) change 



Fig. 5. An example of a Short Characteristic (SC) in 
spherical coordinates (the heavy line connecting U to D 
through P). As in Fig. ||the downstream part of the SC 
(which is only needed to guarantee second order accuracy 
in the formal integral between U and P) is dashed. The 
dotted line shows the complete ray to which the SC be- 
longs. 

along the ray. This means that in addition to spatial in- 
terpolation (see Eq. [I?]) , we are forced to interpolate in /i 
and 4> as well in order to evaluate /„ (U, a;). This is because 
the intensities at points A, B and C are known only for 
a discrete set of directions which are different, in general, 
from (/ii ,(f>i). 

The second, more fundamental difficulty arises be- 
cause in spherical coordinates the concept of upstream 
and downstream boundaries is different from the Carte- 
sian case. Radial infinity is both the upstream and the 
downstream boundary, while in 9 there is no obvious up- 
stream or downstream boundary. If the grid is swept from 
9 = to 9 = 7r or vice versa, one will encounter situ- 
ations in which the intensity at one of the points A, B, 
C is not known before the evaluation the transfer integral 
along the short characteristic (Eq. |l6| ) is performed. An 
example of such a situation is shown in Fig. |^. Interpo- 
lation makes use of the points A, B and C, but since A 
coincides with P, the intensity at the point A = P, has 
not been computed yet. 

4-.1. Extended Short Characteristics 

The problem of unknown intensities can be solved by mod- 
ifying the definition of short characteristics to be the part 
of the ray that connects P, not with just the nearest grid- 
line intersection, but with the nearest R = R^ gridline in- 
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Vx 2 +y 2 

Fig. 6. An example of an Extended Short Characteristic 
(ESC) in spherical coordinates (the heavy line connecting 
U2 to D through P). It is the extended version of the SC 
shown in Fig. ||. The ESC does not have the same problems 
as the SC because for the ESC none of the points A, B or 
C coincides with point P. 

tersection, i.e. the nearest radial shell. Such an "extended 
short characteristic" (ESC) is illustrated in Fig. ^. The 
starting point U of such an ESC will be located either at 
R = Bk-i, B = Bk+i or back at B = Rk- This means 
that in between U and P the ESC may intersect one or 
more grid lines. The point D on the downstream side 
remains the same as for standard SCs. 

By using ESC instead of SC the "problem of un- 
known upstream intensities" can be eliminated. In fact, 
if a proper sweeping scheme is chosen (see Subsection 
4.2), the problem of unknown intensities only occurs in 
those situations when a short characteristic curves back 
onto the same 0-gridline from which it originates, as is 
illustrated in Fig. 0. By extending only those short char- 
acteristics, and leaving the rest truly short, one can also 
avoid unknown intensities in the sweeping scheme. Just 
for notation we call this scheme the Minimally Extended 
Short Characteristics scheme (MESC). MESC is almost 
as accurate as ESC, but ESC is more closely similar to its 
one-dimensional spherical analogues, and is slightly less 
numerically diffusive. 

In the following we denote with D (or £7_i) the sin- 
gle downstream intersection with a grid line, and with 
Ui (i = 1, . . . , m) the multiple intersections upstream of 
P. The point U m is therefore the true upstream starting 
point of the ESC, where the intensity must be found by 
interpolation. In Fig. ^| this is the point t/jj and the ESC 
consists of two segments in this case. 



4-2. The sweeping scheme 

Using the (minimally) extended short characteristics de- 
fined above, we can systematically sweep the grid without 
encountering unknown intensities. We start at the outer 
boundary Boo and integrate inwards only those ESCs for 
which fj,i < 0. The sweeping order in is from = to 
= 7r and then back. 

The intensity at each P = (Bk,Qi) is found for all 
Hi < and (f>j by tracing the ESCs back to their up- 
stream starting point U m . At point U m the values of B, 
0, h and (j) are different from those at P, and will be de- 
noted with R m , m , \x m and <p m . The find the intensity 
at (R m , m , Am, <fim)) by interpolation. We then integrate 
the formal transfer equation along each segment of the 
ESC connecting U m to P, according to Eq. (|l~8|) . This gives 

I V (P) =exp(-T m )I v (U m ) + ^2 exp(-T 4 _i) x 

i=l,m (38) 

where is the depth from P to Ui and the index i denotes 
the quantity evaluated at Ui (e.g. Si = S(Ui); i = refers 
to P and i = — 1 to D). The it's, p's and d's are defined as 
in equations ©, © and (||), but with r replaced by 
(n - Tj_i) and f by (t;_i - Ti_ 2 )- 

The integration is then repeated moving towards 
smaller radii, until the inner boundary at R = R\ is 
reached. Here we can include the contribution of a cen- 
tral source or any other boundary condition. 

Then we start integrating back towards larger radii, 
until we reach the outer edge. By now the radiation field 
on the grid I V (R^, 0; ; Hh 4>j) IS known. 

4-3. Tangent-ray discretization of photon direction 

Equations (|l]) and (|3^) show that fi and change along 
the ESC. If we follow the ESC upstream towards a point 
U m (see Fig. [|) , then the values of these two angles at U rn 
are generally not exactly at the discrete values {fJ-i} and 
{4>j} of the sample of directions. This means that we must 
interpolate not only in space (between A, B and C), but 
also in direction /1, cf>. Although, this is not, in principle, 
a fundamental problem, the use of interpolations should 
be reduced to a minimum to avoid unnecessary numerical 
diffusion. 

Fortunately one can eliminate the interpolation in \x by 
means of a suitable choice of the /i-grid so that all ESCs 
always start and end at one of the {/i^} points. We let the 
//-discretization depend on Rk, 



{Vk,j ; 3 = -TOfc, • • • ,m fc } , 



(39) 



and choose the (ik,j in such a way that for each j = 
—mk, • • • , mfc there is a j' — — m,fe_i, • • • , nik-i such that 



(1 - Mfe,,)^ - (1 - Mfc-i,,')^ 



fe-i 



(40) 
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This choice is based on the fact that the values of \x and 
R along an ESC depend only on each other and on b, as 
can be seen by combining equations (p9h and (|3l|): 



(1 - ,i 2 )i? 2 = b 2 



(41) 



By choosing the Lik,j according to Eq.(flO|), a ray which 
originates at Rk with, say, /i = fik,j and arbitrary <f>, 
reaches any other radius along the path with a value of \x 
which coincides with one of the points of the local \x grid 
there, thus eliminating the need for interpolation. 
A /i-grid that is consistent with Eq.(|40|) is: 



Mfe,±i = ±1 



Rl 



(42) 



in agreement with the angular spacing induced in spherical 
symmetry by the "tangent ray method" (see e.g. Mihalas 



1975; Zane et al. 1996). Actually, it can be easily shown 



that in 1-D spherical symmetry the ESC method, with 
[if. j given by Eq.([l2"|), is fully equivalent to the tangent 
ray method. This is an important feature of the algorithm 
since it then exactly recognizes spherical symmetry. And, 
even in the absence of spherical symmetry, it transports 
radiation outward without any numerical diffusion in R or 
li. 

However, the \x spacing implied by Eq.([42|) has the ten- 
dency to give a poor sampling around \i = 0. This problem 
can be easily solved by introducing some (typically one or 
two) extra points around /z ~ to enhance the angular 
resolution there. Obviously this violates the original pre- 
scription and therefore requires the use of interpolation for 
these extra /z-points, producing a small amount of angular 
diffusion for /i ~ 0. Generally this diffusion is small. 

Unfortunately the interpolation in <f> can never be 
avoided. The (j) angle depends in a complicated way on 
s (see Eq. ^) and it can change rapidly even within one 
element of a ESC. Both first and higher order interpo- 
lation in (j> have been tested in our numerical code. We 
have found that in most cases the <f> diffusion is not very 
large and, in general, influences the solution less than the 
spatial (9) diffusion. 

4-4- Special treatment of radiation near fi ~ 1 

The tangent-ray discretization of [i allows the algorithm 
to accurately conserve radial flux. However, such a choice 
of /i-angles requires a large number of \x points at larger 
radii, typically > k. One cannot make do with a 
smaller number of fx points without facing the risk of loos- 
ing flux. This is illustrated in the following argument. If 
radiation is emitted at a radius i?, an observer at Rk S> R 
can see the radiation from the emitting region even if its 
eyes cannot resolve the source. This is because the ob- 
server's eyes measure the flux and not the intensity. The 
ESC algorithm, on the other hand, deals with intensity, 



and intensity is converted into flux by performing an in- 
tegral over dfidcj). For this integral to be reasonably ac- 
curate, the emitting region must be resolved in and /i, 
leading to the requirement mk ^ k. 

Unfortunately this means that the computational cost 
scales as N\\ if one wishes to extend the span of the ra- 
dial domain. Since the ability to deal with many orders of 
magnitude in radius is crucial to solving transfer problems 
in circumstellar envelopes, this scaling is undesirable. 

An easy way to solve this scaling problem is related to 
the simple observation that all photons with /i ~ 1 follow 
roughly a radially outgoing trajectory and they tend to 
travel more radially (i.e. with \x closer to unity), the fur- 
ther they propagate outwards. In the "radial streaming" 
limit (1 — fx -C 1), Eq.(f40|) becomes approximately 



1 



Hk-lj-l 



^k-l,j-l 



R k-l 

Rl 



(43) 



where ilkj is the solid angle (bounded by fXkj) at Rk- 
Eq.(fl3"l) is just a restatement of the 1/R 2 law which is 
exactly obeyed by a point source and by any radiation in 
the radial streaming limit. 

This property of radially outward radiation makes it 
possible to bundle all /i-points with sufficiently large \x 
into a single collective flux-like bin. The intensity of that 
bin will be treated as the average intensity within that 
collective bin. The idea is to divide the /u-range [—1,1] 
into three parts 



[1, /irs 

inward intensity bin: fj, 

[— fi rs ,fi rs ] intensity samples: fx 

[^rs, 1] outward averaged bin: fx 



-1 
1 . 



This way the number of /z-gridpoints mk at each radius 
can be limited, depending on how close fx rs is to unity. 
We choose a global value for fx rs , and do not allow this to 
differ from one radius to another. 

Because the radial outward bin represents an inte- 
grated intensity, i.e. an average of the true intensity over 
a solid angle il rs = 2n(l — n rs ), it requires a special 
treatment. Let us denote the average intensity in this 
bin as I^(Rk,Qi). The integration formula, Eq.(^8|), for 
/+(i?fc,0;) becomes 



I+(R k ,Q l )=exv(-T) 



R l-l r- 

R 



2 i+(Rk-i,ei) + 



u v S v (Rk-x,®i) + PuS v (Rk,@i) 
d v S v (Rk+i, @i) ■ 



mi,— 1) 



(44) 



This formula simulates the correct 1/R 2 behavior of the 
flux for optically thin media provided fi rs is sufficiently 
close to unity. It reduces to the standard expression when 
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the medium is optically thick. An estimate of the error 
introduced by the assumption of radial streaming can be 
made by comparing equations ( f40| ) and (43), and is of 
order (1 — fi rs )/2. The radiation contained in the n = 1 
bin will not accumulate any numerical interpolation errors 
because it moves strictly along the radial grid lines. 

The inward collective bin will always behave as a real 
intensity, so that the 1/R 2 behavior does not need to be 
simulated. 



based on Long Characteristics. The interpolations used in 
the ESC algorithm introduce numerical diffusion, even in 
the optically thin regime, and this constitutes a potential 
threat to the reliability of the method. In order to test the 
accuracy of the ESC Lambda Operator we have performed 
a series of runs for a number of simple setups, comparing 
the results of the ESC calculation with those obtained by 
means of an exact LC Lambda Operator. Here we present 
the analysis for three such tests. 



4-5. Spectra and images 

Once the iterative part of the transfer has been completed 
and the source function is known, the next step is to pro- 
duce images and spectra. An image is produced by formal 
integration of the source function along long characteris- 
tics through the medium (ray tracing). Each ray repre- 
sents one pixel of the image. One can produce spectra by 
making images at a range of frequencies, and integrating 
these images over the "detector" aperture. 

Here, as in the Lambda iteration, we face resolution 
problems if the source under consideration spans a large 
range in log(i?). The central parts of the image are often 
much brighter than the rest, but cover a much smaller 
fraction of the image. The spectrum may therefore con- 
tain significant contributions of flux from both the cen- 
tral parts and the outer regions of the image. Unless the 
image resolve all spatial scales of the object, the spectra 
produced in such a way are unreliable. 

If a rectangular arrangement of pixels is used, one must 
make sure to use a variable spacing in both x and y, in 
such a way that the small scales around the star are suffi- 
ciently resolved. If one is mainly interested in the images 
themselves, this seems the most reliable and straightfor- 
ward way to go. 

For the production of spectra we propose a different 
approach. Rather than arranging the pixels over a rect- 
angle, we arrange them in concentric rings. The impact 
parameters of the circles are related to the radial grid 
points of the transfer calculation. For a reliable evaluation 
of the spectra it is generally enough to have one circle for 
each Ri, plus some more, about 5, to resolve the central 
region. The number of circles Nf, in each image is there- 
fore roughly the same as the number of radial grid points: 
Ab = Nji + 5. The number of pixels in each circle N v 
is slightly less straightforward to choose, but for reliable 
spectra it is generally sufficient to take N v — 2Nq, where 
Nq is the number of O grid points, counted from pole to 
pole. Using this method, the images automatically resolve 
all relevant scales, while using only a fairly limited number 
of pixels. 

5. Testing the ESC Lambda Operator 

The Extended Short Characteristic implementation of the 
Lambda Operator is not exact, as opposed to the one 



5.1. Optically thick annulus 

The first test problem concerns the determination of the 
radiation field produced by an optically thick, isothermal, 
sharply-edged annuls, bounded by R n < R < Rm and 
Qm < O < 7r — Ojvf. In the actual calculation we have 
taken Rq = 2.86, Rm — 4.83, Oju = 1-26 and the absorp- 
tion a a is given by 

, _ flO 3 inside the annulus , 

^> ) = { O elsewhere. ^ 

For the sake of simplicity all variables are dimensionless 
and the temperature has been taken such that B V (T) = 1. 
We use a spatial grid with 20 radial points, logarithmically 
spaced such that 

Ri+i-Ri = 1 1402; 
Ri 

and the angular grid in 9 consists of 20 equally spaced 
points from pole to pole. Fig. shows the configuration for 
the test problem. Mirror symmetry in the equator reduces 
the number of actual points to 10 in the range < O < 
7r/2. The mesh in photon momentum space consists of 32 
equally spaced points in <fi, covering the range 0-2-7T, and 
41 points in fi, 38 chosen according to Eq. (|40|) , plus fi = 
and fj, = ±0.24 to ensure sufficient resolution at small /i. 
Because of the symmetry, the transfer needs to be solved 
only for the 16 points in the ranges < (f> < ir/2 and 
3tt/2 < 4> < 2tt. 

The radiation field emitted by such a source can be 
semi-analytically determined. As seen by an observer at 
some point P, it is simply the projection of the object on 
the sky of the observer. Since the object is sharply-edged 
and highly optically thick (r ^> 100), its projection will 
be sharply-edged as well. The intensity is simply given by 

, , . _ J for rays missing the object 

v[ ; H, <P) - \B U (T) for rays hitting the object ' 

(46) 

so it is easy to compute the projections on the sky of 
the observer at various positions in space by using some 
independent ray-tracing algorithm or a semi-analytical 
computation of the image. This image can then be com- 
pared to that produced by the ESC and the LC transfer 
algorithms to evaluate the accuracy of the ESC Lambda 
Operator. 
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Fig. 7. The set-up for the test runs. The annulus is shown 
together with the spatial mesh. Here all four quadrants 
are shown, although only the first quadrant needs to be 
computed as a result of mirror symmetry in the equatorial 
plane and cyclindrical symmetry around the polar axis. 
Nine representative grid points, labeled A—H, are marked 
(see Fig. |). 



Let Iesc(^j 4>) an d Ilc{^, 4>) be the intensities, at the 
observer location, computed using the ESC and the LC 
algorithms, for the same discretization in fi and <f). Let, 
moreover, I(fj,, (f>) be the true intensity, which may be 
found by tracing individual rays with very high resolu- 
tion in fi and <j>. We define the standard error of the ESC 
algorithm as 



'ESC 



J (Iesc ~ I) 2 dnd<P 



Jl 2 dfxd<p 

and the error in the mean intensity J as 

J (Iesc - I)dfid(f> _ J ESC - J 



£esc 



J Idfid(f> 



J 



(47) 



(48) 



Similar definitions apply to the LC errors. 

Radiative transfer for this setup has been performed 
using the ESC method. The results are shown in Fig. ^ for 
the 9 gridpoints labeled A — I in Fig. [7| The contours of 
the real images are overplotted. The same calculation has 
been repeated using the LC method. The errors of both 
the ESC and the LC calculations are listed in Table |l|. 
These figures show that the errors of the ESC method are 
not very much greater than those of the LC method (which 
result from the discretization alone) and strcnghtcn the 
reliability of the ESC algorithm. 



Point 


tESC 


a ESC 




ji 


A 


0.060 


0.090 


0.085 


0.116 


B 


0.034 


0.112 


-0.016 


0.134 


C 


0.059 


0.112 


0.001 


0.105 


D 


-0.020 


0.091 


0.002 


0.094 


E 


0.107 


0.096 


-0.032 


0.084 


F 


0.109 


0.103 


0.028 


0.089 


G 


0.029 


0.051 


-0.019 


0.068 


H 


-0.012 


0.056 


-0.046 


0.104 



Table 1. The errors of ESC and the LC algorithms for 
the optically thick test. 



5.2. Optically thin annulus 

Contrary to the optically thick case of the previous sub- 
section, the radiation field from an optically thin annulus 
cannot be determined by a simple analytic formula such 
as Eq.(^). At large radii, however, the mean intensity 
should follow the 1/R 2 law and be independent of 0. We 
can verify if the solution produced by the ESC algorithm 
indeed has this expected behavior. All details are the same 
as in the previous test, with the only difference that now 
the (constant) absorption is taken 10 -6 . 

We focus on the behavior at large radii. In radial 
streaming H(Ri,Qj) ~ J(Ri,Oj) which, for an optically 
thin source, is simply proportional to the volume integral 
of the emissivity j(R, 0) and is independent of 



J(R) 



^ / j(R',@')(R') 2 dR' sine' d& . 
AirR 1 J 



(49) 



Although close to the source the mean intensity is still very 
dependent on 6, at i? > Rm the code should be able to 
recover the correct behaviour J oc R~ 2 at large radii. The 
mean intensity J(Ri, 0j) resulting from the ESC transfer 
calculation is shown in Fig. ^ where it is multiplied with 
R 2 . It can be clearly seen that the while J depends on 
for small radii, it becomes almost independent on as 
the radius increases and that the inverse square law is very 
well reproduced. The dependence of J on at the largest 
radius is shown in Fig. |l0|. Typical errors are < 5%, which 
lies within the error expected from the coarseness of the 
grid. 

5.3. Spherically symmetric test problem 

Although the above tests show that the ESC algorithm 
performs well on its own, they are not enough to prove 
that it will produce accurate and reliable results when 
applied to, for instance, a non-LTE line transfer compu- 
tation. Unfortunately it is not easy to test this, because to 
our knowledge there exists no useful benchmark test case 
yet for 2-D axisymmetric radiative transfer in circumstel- 
lar clouds. 

The least we can do is to test our 2-D algorithm on 
a 1-D spherically symmetric test case, and check the out- 
put against that produced by an independent 1-D transfer 
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Fig. 8. The radiation field I(fj,, cf>) for the test problem of the optically thick annulus at various positions in space (from 
top left to bottom right at the points A ~ H of Fig. In each panel <fi, from to 2tt, is on the cc-axis and /i, from -1 
to 1, on the j/-axis. Tickmarks on both axes are representative of the actual values of /x and 4> used in the calculation. 
The image shows the intensity resulting from the ESC transfer algorithm (the grey scale is such that white corresponds 
to maximum and black to zero). The thin solid curves mark the true contour of the object, computed by ray-tracing 
with high angular resolution. The slightly diffusive nature of the results is a consequence of the interpolations inherent 
to the ESC method. Numerical diffusion remains quite low in almost all cases, with the possible exception of point E. 




Fig. 9. A surface plot of the mean intensty J{R, 6) as pro- 
duced by the ESC algorithm for the optically thin annu- 
lus. In order to show the behavior at large radii, the mean 
intensity is multiplied by R 2 and normalized for conve- 
nience. 



calculation. This way we can at least test two of the spe- 
cial features of the ESC algorithm: the additional /i-points 
close to fj, ~ 0, and the special treatment of the intensity 



1 .2 
1 .0 
8 0.8 

\ 0.6 

o 
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e 

Fig. 10. The computed (solid line) and the exact (dashed 
line) mean intensity for the optically thin annulus, at the 
largest radius as a function of O; errors are < 5%. 



near fi ~ 1. These are features that are not particularly 
related to 2-D, and can therefore also be tested in 1-D just 
as well. 

Our test cloud is a spherically symmetric power law 
model with hydrogen density specified by 



Ro 



(50) 
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where R is the radius in cm, and = 2.0 x 10 7 cm 3 
is the number density at R — Rq = 1.0 x 10 15 cm. We 
take a constant kinetic temperature T^ n (R) = 20 K. The 
abundance of our molecule is also a constant, X mo \(R) = 
iV mol (J?)/iV H2 (Ji) = 1.0 x 10~ 6 . The systematic veloc- 
ity is taken zero. The models are computed in spherical 
coordinates, in the domain i?j n < R < R ou t- We take 
i? in = 1.0 x 10 15 cm and R out = 7.8 x 10 18 cm. At the inner 
boundary we put a reflective boundary condition. The in- 
coming radiation at the outer boundary is the T = 2.728 K 
microwave background radiation. 

We choose a Active 2-level molecule which is specified 

by 



E 2 -Ei = 6.0 cm" 1 = 8.63244 K 

92/ gi = 3.0 
A 21 = 1.0 x 10~ 4 s" 1 
K 2 i = 2.0 x 10~ 10 



S — 1 

cm s 



(51) 
(52) 
(53) 
(54) 



from which the downward collision rate follows: C21 = 
7Vjj 2 _ftT2i s _1 . The total (thermal+turbulent) line width 
atot is a to t = 0.150 kms -1 (see Eq. ||). 

The test problem presented here has high optical depth 
and a very sub-critical density. It is therefore well suited 
to test whether non-LTE effects are properly computed. 

The line transfer is computed in a passband of 40 
frequency points equally spaced between —0.40 kms -1 
and +0.40 kms -1 . The /x angle is discretized using 43 
points, arranged according to the tangent-ray prescrip- 
tion of Eq. ([42|) with 3 additional //-points around /i = 
on each side. Our convergence criterion is simply: 



m&x(Sni/rii) < 1 x 10 



(55) 



at all radii. For the radius we use an equally spaced log- 
arithmic grid with {R l+1 - Ri)/Ri = AR/R = 0.1720. 
We perform 4 runs: Lambda Iteration and Accelerated 
Lambda Iteration with and witout Ng acceleration. 

The results for the upper level population is plot- 
ted in Fig. [ll] and compared to the results obtained 
i ndep endently with SIMLINE written by V. Ossenkopf 
(|1999D . The convergence plots for four different methods 
are shown in Fig. [L2|. 

More test problems, and a more e xtens ive discussion 
of them was presented by Dullemond ( 1999] ). 



6. A simple model for the Egg Nebula 

Now that the algorithm has been tested, we demonstrate 
here how it can be used in practice. Our first example is 
a simple model of the optical appearance of the Cygnus 
Egg Nebula (CRL 2688). This object has been extensively 



studied ever since its discovery by Ney et al. (1975). It is 



a bipolar reflection nebula s urrou nding an F5 supergiant 
of T c ff = 6500K (Crampton 1975 ). At optical wavelength 
it appears as a diffuse bi-lobed nebula with two sharply 
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Fig. 11. The fractional population of the upper level of the 
1-D test problem. The symbols are the solution produced 
by RADICAL, which is our code based on the ESC algo- 
rithm. ALI+Ng were used. The solid line is the solution 
found by the program SIMLINE, which is an independent 
1-D line transfer program written by V. Ossenkopf (1999). 
The difference between the two solutions (normalized to 
the SIMLINE solution) is shown in the lower panels. 



edged "search light b eams" emerging from each of the poles 
(Sahai et al. 1998 ). The lobes are separated by a dark 
equatorial lane which completely obscures the central star. 

The optical emission from this nebula can be un- 
derstood as reflected starlight escap ing fr om the nebul a 
through polar cavities (Latter et al. 1995 , Morris 1981 ). 
It is clear that the "searchlight beams" are due to sin- 
gle scattering of direct starlight by dust grains. The lobes 
are, however, more likely to be the result of multiple 
scattering. 

We will model this multiple-scattering process in 2- 
D with the MESC algorithm, in an attempt to reproduce 
the complex optical appearance of the nebula. Our setup 
consists of an almost spherical wind with a cavity at both 
poles. The density in the cavities is small, but it is still 
high enough to reflect sufficient amounts of starlight. To 
reproduce the twin-beams at both poles, we place a small 
blob of matter at the polar axis in the cavity, causing a 
shadow. A star is placed at the center of the coordinate 
system to illuminate the nebula from within. 

We model only a single frequency in the optical, at 600 
nm. For this reason we refrain from taking actual realistic 
dust opacities, and specify total optical depth and scatter- 
ing albedo instead. The dust density is shown as contour 
plots in Fig. [ll|. The total optical depth at the equator is 
about 60. The ratio of absorption over scattering is 0.3. 
The dust scattering is assumed to be isotropic, which suf- 
fices for the present simplified example. We do not specify 
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Fig. 12. Convergence plot for the 1-D test problem. The 
largest error is plotted against the iteration number for 
four different methods. We define the error as the maxi- 
mum relative difference between the "real" solution and 
the level populations at iteration iV itor . The "real" solu- 
tion was obtained with ALI+Ng and converged down to 
max(8n i /n i ) < 10~ 6 , which is a 100 times stricter con- 
vergence criterion than the last iteration step plotted for 
ALI+Ng in this figure. 
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Fig. 13. The density contours for the dust envelope 
around the central star of the Egg Nebula. Contours are 
logarithmically spaced, separated by factors of \/l0. The 
distance scale on the x- and y-axis is in centimeters. At 
the polar axis one can see a small blob of matter. This 
blob is responsible for the shadow. 



the dust temperature for this setup since thermal emission 
at 600 nm is negligible. 

The simulation was performed by RADICAL, using 
the MESC algorithm, and applying Accelerated Lambda 
Iteration and Ng acceleration. The image subsequently 
produced by formal integration is shown in Fig. The 
model reproduces the searchlight beams and the diffuse 
glow. It also naturally reproduces the intensity difference 
between the north and south lobe. This is a result of the 
slight inclination at which the object is seen. For light 
emerging from the south lobe, the path length through 
the outer regions of the nebula is larger than for the north 
lobe. 

Although the model resemblances the HST image of 
Sahai et al. , it should be noted that the density structure 
that we have used may not be consistent with observations 
at other wavelengths. For example, HCN observations by 
Bicging & Nguyen- Q-Rieu ( 1996 ) seem to rule out the 
presence of a cavity in the wind. Also the rather high 
albedo may be difficult to reconcile with the fact CRL 
2688 is carbon-rich. 



7. A model of line transfer in a collapsing 
protostellar cloud 

In this section it will be demonstrated how the 2-D trans- 
fer algorithm can be used in the observational study of 
low-mass star formation in dense molecular cloud cores. 



Low mass star formation takes place in dense molec- 
ular cloud cor es. A ccording to the spherically symmetric 
model of Shu ( |1977|) , such a core develops a cusp with den- 
sity close to a 1/R 2 powerlaw. Once a gravitational insta- 
bility is trigged, the centeral part of the cloud collapses, 
and forms a star. An expansion wave propagates into the 
cloud towards larger radii, allowing more and more matter 
to fall supersonically down the potential well, and add to 
the protostar's mass. 

Both observational evidence and theoretical argu- 
ments, however, indicate that purely spherical collapse 
is rare. Any slight amount of angular momentum in the 
primordial cloud will cause deviations from sphericity 
as centrifugal forces tend to dominate over radial infall 
deep down the potential well. And even before the col- 
lapse stage these primordi al clo uds often appear to be 
non-spherical (Myers et al. |l99l| ). Theoretical models of 
non-spheric al pro tostellar collapse inclu de, am ong oth- 



ers, Uhich (1976), Cass en fc Moosman ( 1981 ), Terebey 



et al. (|1984 Galli & Shu [1993 



The models of Cassen & Moosman and Uhich (here- 
after CMU) focus on the inner free-falling part of the 
collapsing cloud, and assume that the material originates 
from an originally spherical cloud with some angular mo- 
mentum. Their model is almost spherical at large radii, 
but flattens off closer towards the center, and forms a disk 
near the centrifugal radius. This model was later extended 



by Hartmann et al. (1996, hereafter HCB) to include flat- 



tening of the parent cloud. These models show that the 
inner free-fall part of an initially flattened cloud naturally 
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Fig. 14. The synthetic image of the Egg Nebula, produced by RADICAL. The color coding is reverse logarithmic grey 
scale. The color table was slightly modified to enhance contrast in the lobes, but the modifications remain within 15%. 

tends to form a bipolar cavity, which is often observed 
in YSO. These models are distinctly non-spherical at all 
radii, despite the fact that centrifugal forces only domi- 
nate at small radii. It is therefore evident that fitting such 
models to molecular line observations requires 2-D axi- 
symmetric radiative transfer computations. 

In this section we perform such a calculation, using 
the algorithms of this paper. We solve the non-LTE level 
populations for the first 7 rotational levels of HCO + , and 
compute the predicted single-dish spectra. 



7.1. Description of the model 

In our HCB models we assume that the radius of the ex- 
pansion wave -Race is outside our domain, so we shall con- 
fine our study to the free-fall inner region of the collapsing 
sheet-like molecular cloud. We assume that matter in the 
parent cloud had a small amount of rotation in the plane 



of the sheet before it collapsed. According to the HCB 
model, the ve locity field of the gas is given by the formu- 
lae of Ulrich fll976|) : 



vr 



VQ 




(56) 
(57) 

(58) 



where [i = cos(0) and = cos(0o)- The angle 6o is the 
O-coordinate that the gas parcel had when it started its 
free-fall at large radius. For a given (R, 0) the value of 6o 
can be found from 



R 



1 



JL 

Mo 



(1 " Mo) , 



(59) 
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Fig. 16. A zoom-in of the density structure of the 77 = 3 
case of the HCB model. The figure on the right is a zoom- 
in to the scale of the centrifugal radius. The bar represents 
the disk. 

where R c is the centrifugal radius, i.e. the radius at which 
centrifugal forces equal gravity. This is the outer radius of 
the disk that is formed as a result of the rotation. 

The density of the gas for the HCB model is given by 



Phcb(R, ©) 



M 



"47r(GMi? 3 ) 1 /2 tanh(77) 



1 



Mo 



-1/2 



Mo 



2^R C 
R 



_ x (60) 



where n is a dimcnsionlcss flattening parameter, which is 
roughly equal to the ratio of the accretion radius R acc to 
the sheet thickness H . HCB argue that this value must be 
somewhere in between 77 = and 77 = 77 max — 4. For rj = 
the CMU models are reproduced. Density contours of this 
free-falling envelope, for different flattening parameters, 
are shown in Fig. [l5| The centrifugal radius of the infalling 
envelope is at R c = 100 AU. We place a thin disk with a 
radius of i?disk — R c = 100 AU at the equator. A zoom-in 
of the density distribution, down to the scale of the disk, 
is shown in Fig. |l6|. We will ignore any emission from 
the disk, and merely treat it as a light-blocking boundary 
condition at the equator. 

7.2. Non-LTE line transfer 

We compute the line transfer problem for four HCB mod- 
els, with flattening parameter 77 = 0,1,2,3 for models 
1,2,3,4 respectively. The adopted valued for the accre- 
tion rate and the turbulent line width are M — 2.4 x 
lO~ 6 M yr _1 and a tur b = 0.25 kms -1 . 

For these models we compute the non-LTE line trans- 
fer problem of the lowest few rotational levels of HCO + , 
including the effects of the moving medium. We assume a 
cosmic background (CMB) continuum as incident radia- 
tion at the outer edge of the computational domain. Dust 
emission and opacity are neglected in the line transfer, 
which is justified for the lower-lying HCO + lines because 
the nebula is optically thin to dust in the millimeter and 
sub- millimeter, and radiative pumping by dust continuum 



is not important for HCO + . Also, we need not include the 
dust emissivity in the computation of the emerging spec- 
tra, since we shall show only the spectra with the dust- 
and CMB-continuum removed. The radiative transfer is 
computed within a range of R € [70, 10 4 ] AU. As an inner 
boundary we have a vacuum. The cross sections for H2- 
HCO + collisional transitions were taken from Monteiro 
( |1985D and Green ( |1975[ ). We adopt an HCO+ abundance 
of 2 x 10 -9 The gas temperature is taken to be T = 20K 
throughout the cloud. 

We perform the non-LTE line transfer for all four mod- 
els. The resulting non-LTE level populations for model 3 
are shown in Fig. [l?], and the corresponding excitation 
temperatures are shown in Fig. [l^. One can see that at the 
equator (0 = 7r/2) the levels are almost thermalized, ex- 
cept at large radii. This is due to the much larger density 
at the equator than at the pole. The drop in excitation 
temperature at large radii is a result of the decoupling 
of radiation and matter. At large radii the level popu- 
lations will be strongly influenced by the cosmic back- 
ground radiation. Another interesting phenomenon occurs 
at small radii near the pole: the excitation temperature 
exceeds the gas temperature. This effect was discussed 
by Leung & Liszt (1976) for the CO 1-0 transition. It 



can be understood as resulting from an overpopulation of 
the J = 1 level due to the large ratio of radiative rates 
(A 21 /A w ~ 10). 

Once the level populations have been computed, the 
line spectra are produced. The spectra are centered on 
the origin of the object. First the circular images are pro- 
duces in a range of frequencies. This circular rendering of 
the images ensures that no details at large or small radii 
are missed, and thus that no flux is accidently lost. The 
antenna temperatures are then computed by integrating 
the images, after they have been multiplied by the beam 
pattern centered on the origin of the object. We use an 
"Airy" beam, with a beam size corresponding to a single 
dish of 15 m diameter. The object is placed at 140 par- 
sec distance. The spectra of the four models, computed 
for the first four radiative transitions at three different 
inclinations, are shown in Fig. [l9| 

From the spectra one can clearly see the effect of flat- 
tening of the HCB cloud, in particular for models 3 and 4. 
At near pole-on inclination (5°) hardly any self-absorption 
is seen in these models, because one looks straight into 
the "cavity". At near edge-on inclination (85°) the "torus" 
blocks the central regions from view near line center, which 
results in the clear self-absorption features seen in the line 
shapes. A similar manifestation of the non-spherical sym- 
metry of the cloud has been discussed recently by Van der 
Tak et al. 199£. An interesting feature of the line spec- 



tra of models 3 and 4 is that the edge-on line profiles are 
wider than the pole-on profiles. This can be attributed to 
the fact that the density and the excitation temperature 
is lower at the pole than at the equator. At pole-on incli- 
nation the high density equatorial matter will emit near 
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Fig. 15. The density structure of the flattened collapsing cloud model of HCB, for varying flattening parameter r\. 
Scales are in AU. The centrifugal radius is at R c = 100 AU. The effect of rotation can be seen most clearly in the 
rj = case, which is spherical at large radius, but flattens off near R = R c . 




Fig. 17. The level populations for model 3 for the first 
four rotational levels of HCO + as a function of and R. 
Note that the centrifugal radius (which is the disk outer 
edge) is at R c = 1.5 x 10 15 cm. 

line-center instead of in the line wings, thus making the 
line profile narrower. 

The asymmetry between the red-shifted and the blue- 
shifted peaks are typical for protostellar collapse. The ro- 
tation is hardly seen in these spectra. This is because the 
rotational velocity is everywhere much smaller than the 
free-fall velocity, except at very small radii where the emis- 
sion barely contributes to the single-dish spectra shown 
here. 

8. Conclusions 

Numerical radiative transfer modeling on desktop work- 
stations is extremely cheap. Not doing so in cases where 
this is possible would mean an enormous waste of valuable 
information that lies encoded in observed data. However 
the success of such modeling depends on the algorithms 



Fig. 18. The excitation temperatures of the lowest four 
rotational transitions of HCO + for model 3, as a func- 
tion of and R. They were deduced from the level pop- 
ulations shown in Fig. using the formula T C x(r/) = 
{hvij/k)/log(njgi/nigj). 

that are available. We have developed a robust and accu- 
rate method, called the "extended short characteristics" 
(ESC) method, by which complicated 2-D axi-symmetric 
multi-frequency radiative transfer calculations can be per- 
formed. By using spherical coordinates, this method can 
accurately treat circumstellar envelopes and disks from 
the stellar surface all the way up to parsec scale, without 
the need of grid refinement. By making a special choice 
of discrete photon angles and bundling 'almost-radially- 
moving' rays into a single bin, the conservation of radial 
flux can be guaranteed even over many orders of magni- 
tude in radius, and without excessive computational cost. 

The ESC method, and a slight variation called MESC, 
forms the core of a multi-purpose 2-D radiative transfer 
code called RADICAL. We have tested the ESC/MESC 
algorithm on a simple test problem which we described 
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Fig. 19. The single-dish HCO+ line 
5° (near pole-on), dashed line is 45° 
object is located at 140 pc. 



spectra of the HCB models, shown at three different inclinations: dotted line is 
and solid line is 85° (near edge-on). A dish diameter of 15 m was taken, and the 



in this paper. We have also verified that the 2-D algo- 
rithm, when applied to a 1-D spherically symmetric prob- 
lem, indeed reproduces what an independent 1-D algo- 
rithm would produce for the same problem. The errors 
remained within a few percent, the exact value of which 
depends on the grid resolution. 

The ESC/MESC algorithm is designed for a variety 
of applications. We have demonstrated in this paper how 
the method can be used for the problem of dust scatter- 
ing in a bipolar proto-planetary nebula and the problem 
of non-LTE line transfer in a collapsing cloud. But the 
method can also easily be applied to other radiative pro- 
cesses, such as dust continuum emission with radiative 
equilibrium for the dust grains, thermal Bremsstrahlung, 
electron scattering and even Comptonization in hot plas- 



mas. The ESC/MESC algorithms are therefore a useful 
tool for a wide range of astrophysical problems. 
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